function feature = colorPyramid(grayLayer, ALayer, BLayer, numBins,numLevels)
pyramid = cell(2^(numLevels-1),2^(numLevels-1),numLevels);

[sx sy]=size(grayLayer); % assume the three layers are of the same size;
% for level = numLevels

% for level numLevels
for b_x = 1:2^(numLevels-1)
    for b_y = 1:2^(numLevels-1)
        x_block_size = ceil(sx*0.5^(numLevels-1));
        y_block_size = ceil(sy*0.5^(numLevels-1));
        x = x_block_size*(b_x-1)+1:min(sx,x_block_size*b_x);
        y = y_block_size*(b_y-1)+1:min(sy,y_block_size*b_y);
        c_gray = imhist(grayLayer(x,y),numBins(1));
        c_A = imhist(ALayer(x,y),numBins(2));
        c_B = imhist(BLayer(x,y),numBins(3));
        counts = [c_gray/sum(c_gray); c_A/sum(c_A); c_B/sum(c_B)]; 
	counts = max(counts,1e-6);
        pyramid{b_x,b_y,numLevels}=counts;
    end
end

    
for iLevel = numLevels-1:-1:1
    for bx = 1:2^(iLevel-1)
        for by = 1:2^(iLevel-1)
            counts = pyramid{bx*2-1,by*2-1,iLevel+1};
            counts = counts + pyramid{bx*2,by*2-1,iLevel+1};
            counts = counts + pyramid{bx*2-1,by*2,iLevel+1};
            counts = counts + pyramid{bx*2,by*2,iLevel+1};
            pyramid{bx,by,iLevel}=counts/4;
        end
    end
end

feature = [];
for iLevel = numLevels:-1:1
        for b_x = 1:2^(iLevel-1)
        for b_y = 1:2^(iLevel-1)
            feature = [feature;pyramid{b_x,b_y,iLevel}];
        end
        end
end

if numLevels>=3
    outter_hist = zeros(sum(numBins),1);
    inner_hist = outter_hist;
    for j = 1:4
	outter_hist = outter_hist + pyramid{1,j,3} + pyramid{4,j,3};
    end
    outter_hist = outter_hist + pyramid{2,1,3} + pyramid{2,4,3};
    outter_hist = outter_hist + pyramid{3,1,3} + pyramid{3,4,3};
    outter_hist = outter_hist/12;
    outter_hist = outter_hist(1:numBins(1));

    corner_hist = pyramid{1,1,3}+pyramid{4,4,3} + pyramid{1,4,3} + pyramid{4,1,3};  
    corner_hist = corner_hist(1:numBins(1))/4; 


    inner_hist = inner_hist + pyramid{2,2,3} + pyramid{2,3,3} + pyramid{3,2,3} + pyramid{3,3,3};
    inner_hist = inner_hist/4;
    inner_hist = inner_hist(1:numBins(1));

    feature = [feature; reshape(min(inner_hist,corner_hist),[],1)];
end
